#################################################################
## Tyler Pratt ##################################################
## File to Replicate Statistical Results  for ################### 
## "Deference and Hierarchy in International Regime Complexes" ##
## June 22, 2017 ################################################
#################################################################
setwd()

## Read in data 
#Def <- read.csv("~/Dropbox/IO SNA/IO_Submission/Deference_Replication.csv")
Def <- read.csv("Deference_Replication.csv")

#########################
## Descriptive Figures ##
#########################

## Figure 1: Institutional Deference with UN General Assembly, 2010
UN.2010 <- Def[Def$regime=="CT" & Def$IO2 == "UN"  & Def$year == 2010,]
UN.coords <- UN.2010$deference_prop
names(UN.coords) <- as.character(UN.2010$IO1)

UN.2010.o <- Def[Def$regime=="CT" & Def$IO1 == "UN"  & Def$year == 2010,]
UN.coords.o <- UN.2010.o$deference_prop
names(UN.coords.o) <- as.character(UN.2010.o$IO2)

UN.coords <- UN.coords[intersect(names(UN.coords), names(UN.coords.o))]
UN.coords.o <- UN.coords.o[intersect(names(UN.coords), names(UN.coords.o))]

mat <- rbind(UN.coords, UN.coords.o)
mat <- mat[,order(mat[1,], decreasing=T)]
mat2 <- mat[,c("ICAO", "ARF", "NATO", "OAS", "UNSC",  "EU", "FATF", "G8", "SAARC",
                "SCO", "COE", "OSCE")]

barplot(mat2, beside = TRUE,
        ylab = "Deference (% of Policy Documents)", ylim = c(0, 1),
        legend = c("Deference to UN General Assembly", 
                   "Deference from UN General Assembly"))


## Figure 2: Networks of Institutional Deference
library("igraph")

# Intellectual Property Rights network
Def.IPR <- Def[Def$regime=="IPR",]
orgs <- intersect(as.character(Def.IPR$IO1), as.character(Def.IPR$IO2))
IPR <- matrix(nrow=length(orgs), ncol=length(orgs))
rownames(IPR) <- colnames(IPR) <-orgs
for(i in 1:nrow(IPR)){
  sub <- Def.IPR[Def.IPR$IO1 == rownames(IPR)[i],]
  avg.def <- tapply(sub$deference_prop, sub$IO2, mean, na.rm=T)
  IPR[i,] <- avg.def[match(rownames(IPR), names(avg.def))]
}
diag(IPR) <- 0
IPR.G <- graph.adjacency(IPR, weighted = TRUE, mode="directed")

set.seed(21112)
layout.IPR <- layout.auto(IPR.G)
layout.IPR[1,] <- c(6.45, -3.1)
layout.IPR[2,] <- c(6.904524, -2.99435)
layout.IPR[5,] <- c(6.95, -4.3)
layout.IPR[4,] <- c(6.423321, -4.390889)
layout.IPR[6,] <- c(6.349443, -3.806356)

plot.igraph(IPR.G, edge.width = (E(IPR.G)$weight+0.1)*12,
            edge.color = "gray47",
            vertex.frame.color = "dark grey",
            vertex.color = "light green",
            vertex.label = V(IPR.G)$name,         
            vertex.label.color = "black",
            vertex.label.cex = 1.7,
            vertex.size = (colSums(IPR)+.8)*30,
            edge.arrow.size = 1.6, edge.curved = T,
            layout = layout.IPR)

## Election Monitoring network
Def.EM <- Def[Def$regime=="EM",]
orgs <- intersect(as.character(Def.EM$IO1), as.character(Def.EM$IO2))
EM <- matrix(nrow=length(orgs), ncol=length(orgs))
rownames(EM) <- colnames(EM) <-orgs
for(i in 1:nrow(EM)){
  sub <- Def.EM[Def.EM$IO1 == rownames(EM)[i],]
  avg.def <- tapply(sub$deference_prop, sub$IO2, mean, na.rm=T)
  EM[i,] <- avg.def[match(rownames(EM), names(avg.def))]
}
diag(EM) <- 0
EM.G <- graph.adjacency(EM, weighted = TRUE, mode="directed")

set.seed(21112)
layout.EM <- layout.auto(EM.G)
layout.EM[1,] <- c(6.51385, 4.097679)
layout.EM[4,] <- c(7.365623, 4.293747)
layout.EM[9,] <- c(7, 3.698653)
layout.EM[5,] <- c(6.516630, 4.791457)
layout.EM[2,] <- c(4, 3.774704)
layout.EM[8,] <- c(4.7, 4.1)

plot.igraph(EM.G, edge.width = (E(EM.G)$weight+0.1)*12,
            edge.color = "gray47",
            vertex.frame.color = "dark grey",
            vertex.color = "light blue",
            vertex.label = V(EM.G)$name,         
            vertex.label.color = "black",
            vertex.label.cex = 1.7,
            vertex.size = (colSums(EM)+.75)*25,
            edge.arrow.size = 1.6, edge.curved = T,
            layout = layout.EM)

## Counterterrorism network
Def.CT <- Def[Def$regime=="CT",]
orgs <- intersect(as.character(Def.CT$IO1), as.character(Def.CT$IO2))
CT <- matrix(nrow=length(orgs), ncol=length(orgs))
rownames(CT) <- colnames(CT) <-orgs
for(i in 1:nrow(CT)){
  sub <- Def.CT[Def.CT$IO1 == rownames(CT)[i],]
  avg.def <- tapply(sub$deference_prop, sub$IO2, mean, na.rm=T)
  CT[i,] <- avg.def[match(rownames(CT), names(avg.def))]
}
diag(CT) <- 0
CT.G <- graph.adjacency(CT, weighted = TRUE, mode="directed")

set.seed(21112)
layout.CT <- layout.auto(CT.G)
layout.CT[16,] <- c(5.1, 3.9)
layout.CT[15,] <- c(3.45, 4.75)
layout.CT[14,] <- c(4.3, 5.6)
layout.CT[13,] <- c(5.65, 5.2)
layout.CT[12,] <- c(1.4, 4.6)
layout.CT[6,] <- c(3.1, 3.8)
layout.CT[2,] <- c(3.6, 5.64)
layout.CT[9,] <- c(2, 4.1)
layout.CT[4,] <- c(1.1, 4.018138)
layout.CT[10,] <- c(4.9, 5.471045)
layout.CT[7,] <- c(1.9, 5.25)
layout.CT[5,] <- c(1.7, 3.5)
layout.CT[11,] <- c(2.4, 3.3)
layout.CT[1,] <- c(3.64324, 3.156588)
layout.CT[8,] <- c(3.8, 3.45)
layout.CT[3,] <- c(6.15, 4.75)
layout.CT[1,] <- c(3.55, 3.2)

plot.igraph(CT.G, edge.width = (E(CT.G)$weight+0.1)*9,
            edge.color = "gray47",
            vertex.frame.color = "dark grey",
            vertex.label = V(CT.G)$name,         
            vertex.label.color = "black",
            vertex.label.cex = 1.5,
            vertex.size = (colSums(CT, na.rm=T)+1)*10,
            edge.arrow.size = 1.5, edge.curved = T,
            layout = layout.CT)


## Figure 3: Division of Labor with OSCE, 2008
# read in data from fitted counterterorism IO topic model (see Deference_TopicModel.R for estimation code)
subissues <- read.csv("CT.subissues.csv")
subissues$aml <- subissues$aml1 + subissues$aml2 + subissues$aml3
subissues$other <- subissues$coop + subissues$dip1 + subissues$dip2 + subissues$dip3_aml

rows1 <- which(subissues$IGO %in% c("OSCE", "SCO") & subissues$year == "2012")
mat1 <- t(subissues[rows1, c("other", "hr", "nuclear_trans", "aml", "crim")])
rows2 <- c(which(subissues$IGO == "OSCE" & subissues$year == "2012"),
           which(subissues$IGO == "EU" & subissues$year == "2012"))
mat2 <- t(subissues[rows2, c("other", "hr", "nuclear_trans", "aml", "crim")])

library(RColorBrewer)
sequential <- brewer.pal(5, "YlGnBu")

barplot(cbind(mat1, rep(0, 5), mat2), ylab = "Topic Distribution",
        names.arg = c("OSCE", "SCO", "", "OSCE", "EU"),
        col = sequential, xlim = c(0, 8))
legend("topright", legend = c("Criminal Justice", "Terrorism Finance", 
                              "Transport & Nuclear", "Human Rights", "Other"), 
       cex = 0.9, fill = sequential[5:1], title = "Topics")

## Compare DOL score for the two IO pairs in 2012
sum(abs(subissues[rows1[1], 4:13] - subissues[rows1[2], 4:13])) # 0.297 DOL score for OSCE, SCO
sum(abs(subissues[rows2[1], 4:13] - subissues[rows2[2], 4:13])) # 1.170 DOL score for OSCE, EU

## Calculate DOL scores for all IO pairs in CT regime complex
DOL <- as.data.frame(matrix(nrow=0, ncol=5))
colnames(DOL) <- c("IO1", "IO2", "year", "DOL", "Def_ind")
for(i in 1999:2013){
  IOs <- unique(subissues$IGO[subissues$year==i])
  if(length(IOs)>1){
    IO.pairs <- t(combn(IOs, 2))
    subissues.yr <- subissues[subissues$year==i,]
    DOL.yr <- cbind(as.character(IO.pairs[,1]), as.character(IO.pairs[,2]),
                    rep(i, nrow(IO.pairs)), rep(NA, nrow(IO.pairs)), rep(NA, nrow(IO.pairs)))
    for(j in 1:nrow(IO.pairs)){
      DOL.score <- sum(abs(subissues.yr[subissues.yr$IGO==as.character(IO.pairs[j,1]),4:13] - 
                           subissues.yr[subissues.yr$IGO==as.character(IO.pairs[j,2]),4:13]))
      DOL.yr[j,4] <- DOL.score
      DOL.yr[j,5] <- ifelse(any(Def.CT$deference_prop[Def.CT$IO1 %in% c(DOL.yr[j,1:2]) & 
                                                     Def.CT$IO2 %in% c(DOL.yr[j,1:2])] > 0), 1, 0)
    }
    colnames(DOL.yr) <- colnames(DOL)
    DOL <- rbind(DOL, DOL.yr)
  }
}
DOL$DOL <- as.numeric(paste(DOL$DOL))
mean(DOL$DOL) # avg. DOL score of 1.00
t.test(DOL$DOL[DOL$Def_ind==1], DOL$DOL[DOL$Def_ind==0]) # p < 0.001





#####################
## Empirical Tests ##
#####################

## Function for clustering SEs
cl   <- function(dat,fm, cluster){ 
  require(sandwich, quietly = TRUE)
  require(lmtest, quietly = TRUE)
  M <- length(unique(cluster))
  N <- length(cluster)
  K <- fm$rank
  dfc <- (M/(M-1))*((N-1)/(N-K))
  uj  <- apply(estfun(fm),2, function(x) tapply(x, cluster, sum));
  vcovCL <- dfc*sandwich(fm, meat=crossprod(uj)/N)
  coeftest(fm, vcovCL) }


## Table 1: Deference and Divison of Labor

Def.CT <- Def[Def$regime=="CT",]

DOL_CT_bl <- lm(DOL_change ~ deference_prop_lag + BindingIO + 
               Mem.Overlap + abs(majors_diff) + UNIP.diff.avg, data = Def.CT)
dat <- Def.CT[-DOL_CT_bl$na.action,]
cl(dat, DOL_CT_bl, as.character(dat$Dir.DyadID)) # deference effect in CT: 0.039*

DOL_CT2 <- lm(DOL_change ~ deference_prop_lag + BindingIO + 
               Mem.Overlap + abs(majors_diff) + UNIP.diff.avg + 
               as.factor(DyadID), data = Def.CT)
dat <- Def.CT[-DOL_CT2$na.action,]
cl(dat, DOL_CT2, as.character(dat$Dir.DyadID)) # deference effect in CT: 0.071**

DOL_all <- lm(DOL_change ~ deference_prop_lag + BindingIO + 
                Mem.Overlap + abs(majors_diff) + UNIP.diff.avg + 
                as.factor(DyadID), data = Def)
dat <- Def[-DOL_all$na.action,]
cl(dat, DOL_all, as.character(dat$Dir.DyadID)) # deference effect across all: 0.047**


## Table 2: Determinants of Deference

Det.bl <- glm(deference_prop ~ 
                majors_diff + IO2_weight + 
                IO2_Technical + Binding_Technical + DecisionMaking_Ease_Dir +
                regime + t + I(t^2) + I(t^3),
              data = Def, family = binomial(link = "logit"))
cl(Def, Det.bl, as.character(Def$Dir.DyadID))

Det.1 <- glm(deference_prop ~ 
               majors_diff + IO2_weight + 
               IO2_Technical + Binding_Technical + DecisionMaking_Ease_Dir + 
               Mem.Overlap + Nested + 
               IO2_funds_log + 
               UNIP.diff.avg + 
               regime + 
               regional_global + 
               prev_def_otherdirection + 
               t + I(t^2) + I(t^3),
             data = Def, family = binomial(link = "logit"))
dat <- Def[-Det.1$na.action,]
cl(dat, Det.1, as.character(dat$Dir.DyadID))

# This model generates warning: fitted probabilities numerically 0 or 1, so I use bayesian implementation below
Det.2 <- glm(deference_prop ~  
               majors_diff + IO2_weight + 
               IO2_Technical + Binding_Technical + DecisionMaking_Ease_Dir + 
               Mem.Overlap + Nested + 
               IO2_funds_log + 
               UNIP.diff.avg + 
               IO2_Technical + 
               regional_global + 
               prev_def_otherdirection + 
               regime + 
               as.factor(DyadID) + 
               t + I(t^2) + I(t^3),
             data = Def, family = binomial(link = "logit"))
dat <- Def[-Det.2$na.action,]
cl(dat, Det.2, as.character(dat$Dir.DyadID))

library(arm)
Det.2b <- bayesglm(deference_prop ~ 
                     majors_diff + IO2_weight + 
                     IO2_Technical + Binding_Technical + DecisionMaking_Ease_Dir + 
                     Mem.Overlap + Nested + 
                     IO2_funds_log + 
                     UNIP.diff.avg + 
                     IO2_Technical + 
                     regional_global + 
                     prev_def_otherdirection + 
                     regime + 
                     as.factor(DyadID) + 
                     t + I(t^2) + I(t^3),
                   data = Def, family = binomial(link = "logit"))
dat <- Def[-Det.2b$na.action,]
cl(dat, Det.2b, as.character(dat$Dir.DyadID))


## Import bootstrapped results for model 1 to estimate substantive effects
# (see bootstrap_deference.R)
boot.results <- read.csv("bootBL_Def2017.csv", row.names=1)

## Figure 4: Substantive Effects of Power, Efficiency Variables on Institutional Deference
plot(colMeans(boot.results)[c(5,4,3,2,1)], 1:5, bty="n", pty="sq",
     yaxt = "n", ylim = c(0.4, 5.6), ylab="",
     xlab = "Change in Deference Intensity", pch = 16, cex = 1,
     xlim = c(-0.045, 0.11),
     col = c(rep("blue", 3), rep("red", 2)))
abline(v=0, lty=2)
lines(x=quantile(boot.results$Decision, probs = c(0.025, 0.975)), y=c(1,1), col = "blue", lwd=2)
lines(x=quantile(boot.results$IO2_weight, probs = c(0.025, 0.975)), y=c(2,2), col="blue", lwd=2)
lines(x=quantile(boot.results$IO2_tech, probs = c(0.025, 0.975)), y=c(3,3), col = "blue", lwd=2)
lines(x=quantile(boot.results$IO2_weight, probs = c(0.025, 0.975)), y=c(4,4), col="red", lwd=2)
lines(x=quantile(boot.results$majors_diff, probs = c(0.025, 0.975)), y=c(5,5), col = "red", lwd=2)
text(x=colMeans(boot.results)[5], y=1.15, "Decision-Making Difference", col = "blue", cex=0.9)
text(x=colMeans(boot.results)[4], y=2.15, "Binding-Technical Pair", col = "blue", cex=0.9)
text(x=colMeans(boot.results)[3], y=3.15, "Technical IO", col = "blue", cex=0.9)
text(x=colMeans(boot.results)[2], y=4.15, "Weighted Voting", col = "red", cex=0.9)
text(x=colMeans(boot.results)[1], y=5.15, "Major Power Difference", col = "red", cex=0.9)


## Robustness

## linear model
OLS.bl <- lm(deference_prop ~ 
                majors_diff + IO2_weight + 
                IO2_Technical + Binding_Technical + DecisionMaking_Ease_Dir +
                regime + t + I(t^2) + I(t^3),
              data = Def)
cl(Def, OLS.bl, as.character(Def$Dir.DyadID))

## probit
probit.bl <- glm(deference_prop ~ 
               majors_diff + IO2_weight + 
               IO2_Technical + Binding_Technical + DecisionMaking_Ease_Dir +
               regime + t + I(t^2) + I(t^3),
             data = Def, family = binomial(link="probit"))
cl(Def, probit.bl, as.character(Def$Dir.DyadID))

## binomial
bin.bl <- glm(cbind(deference_successes, deference_failures) ~ 
                majors_diff + IO2_weight + 
                IO2_Technical + Binding_Technical + DecisionMaking_Ease_Dir +
                regime + t + I(t^2) + I(t^3),
              data = Def, family = binomial)
cl(Def, bin.bl, as.character(Def$Dir.DyadID))

## full scale of outcome variable
scale.bl <- lm(coord.avg ~ 
                  majors_diff + IO2_weight + 
                  IO2_Technical + Binding_Technical + DecisionMaking_Ease_Dir +
                  regime + t + I(t^2) + I(t^3),
                data = Def)
cl(Def, scale.bl, as.character(Def$Dir.DyadID))

## GDP in lieu of major power
GDP.bl <- glm(deference_prop ~ 
                GDP.Diff.l + IO2_weight + 
                IO2_Technical + Binding_Technical + DecisionMaking_Ease_Dir +
                regime + t + I(t^2) + I(t^3),
              data = Def, family = binomial(link="logit"))
cl(Def, GDP.bl, as.character(Def$Dir.DyadID))

## NMC in lieu of major power
NMC.bl <- glm(deference_prop ~ 
                CINC.diff + IO2_weight + 
                 IO2_Technical + Binding_Technical + DecisionMaking_Ease_Dir +
                 regime + t + I(t^2) + I(t^3),
               data = Def, family = binomial(link="logit"))
dat <- Def[-NMC.bl$na.action,]
cl(dat, NMC.bl, as.character(dat$Dir.DyadID))


## remove instituitonal design variables
nodesign.bl <- glm(deference_prop ~ 
                     majors_diff + IO2_weight + 
                    regime + t + I(t^2) + I(t^3),
              data = Def, family = binomial(link="logit"))
cl(Def, nodesign.bl, as.character(Def$Dir.DyadID))




